more-predictions

Cardiac Output

Cardiac output is the product of heart rate (hr) and stroke volume (sv).

We will use a parametric bootstrap to make predictions of cardiac output as a function of log10(body mass) and habitat based on our random-effect models for hr and sv.

n_boot <- 1000

We will make sets of many (1000) predictions for each of the species that occur in both the hr and sv data sets.

# note NAMES HAVE TO BE SAME as passed to glmmTMB() when models were fitted (geesh)
mammal_hr <- readRDS('fitted-models/hr-data.RDS')
hr_mod <- readRDS('fitted-models/hr-re-model.RDS')
mammal_sv <- readRDS('fitted-models/sv-data.RDS')
sv_mod <- readRDS('fitted-models/sv-re-model.RDS')

cardiac_species <- intersect(pull(mammal_hr, animal), pull(mammal_sv, animal))

The two datasets have 24 species in common.

Make dataset for which to make predictions, containing these species.

cardiac_hr <- mammal_hr %>%
  dplyr::filter(animal %in% cardiac_species) %>%
  dplyr::select(animal, order, genus, species, log10.body.mass, log10hr, habitat) %>%
  mutate(across(where(is.factor), factor))

cardiac_sv <- mammal_sv %>%
  dplyr::filter(animal %in% cardiac_species) %>%
  dplyr::select(animal, order, genus, species, log10.body.mass, log10sv, habitat)

For each row of the dataset, make a set of 1000 model predictions of cardiac output. These predictions will include uncertainty in the model parameters (intercept and slopes) as well as variation attributed to order/genus/species. (Residual variation is not included). The mass estimates are different in the different datasets; here we use the mass estimate from the dataset to which the model (hr or sv) was originally fitted.

First we define custom functions to make the predictions

predict_hr <- function(.){
  predict(.,
          newdata = cardiac_hr,
          re.form = NULL,
          type = 'response')
}

predict_sv <- function(.){
  predict(.,
          newdata = cardiac_sv,
          re.form = NULL,
          type = 'response')
}

Do the bootstrap(s)

if (!file.exists('fitted-models/hr_boot.RDS')){
hr_boot <- bootMer(hr_mod, 
                   FUN = predict_hr,
                   nsim = n_boot,
                   re.form = NULL,
                   type = 'parametric'
                   )

saveRDS(hr_boot, 'fitted-models/hr_boot.RDS')
}else{
  hr_boot <- readRDS('fitted-models/hr_boot.RDS')
}

if (!file.exists('fitted-models/sv_boot.RDS')){
sv_boot <- bootMer(sv_mod, 
                   FUN = predict_sv,
                   nsim = n_boot,
                   re.form = NULL,
                   type = 'parametric')
saveRDS(sv_boot, 'fitted-models/sv_boot.RDS')
}else{
  sv_boot <- readRDS('fitted-models/sv_boot.RDS')
}

Add boot results to datasets, combine to one dataset, and compute predicted cardiac output (and its median and 95% bootstrap CI):

cardiac_hr <- cardiac_hr %>%
  mutate(boot_pred_hr = unlist(apply(t(hr_boot$t), 1, list), recursive = FALSE))

cardiac_sv <- cardiac_sv %>%
  mutate(boot_pred_sv = unlist(apply(t(sv_boot$t), 1, list), recursive = FALSE))
  
cardiac_preds <- inner_join(cardiac_hr, cardiac_sv,
                         by = c('animal', 'order', 'genus', 'species', 'habitat')) %>%
  mutate(boot_pred_cardiac = map2(boot_pred_hr, boot_pred_sv, function(.x, .y) mapply(prod, 10^.x, 10^.y))) %>%
  mutate(boot_median_cardiac = map_dbl(boot_pred_cardiac, median, na.rm = TRUE),
         boot_lo = map_dbl(boot_pred_cardiac, quantile, probs = 0.025, na.rm = TRUE),
         boot_hi = map_dbl(boot_pred_cardiac, quantile, probs = 0.975, na.rm = TRUE),
         mean_mass = ((10^log10.body.mass.x + 10^log10.body.mass.y) / 2) )
# units: ml/stroke * beats/min

plot results

cardiac_mass <- gf_point(boot_median_cardiac ~ mean_mass, data = cardiac_preds,
         color = ~habitat,
         text = ~animal) %>%
  gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = cardiac_preds,
              color = ~habitat) %>%
  gf_labs(title = 'Expected Cardiac Output\nfor the species modeled',
          y = 'Predicted Cardiac Output (mL/min)',
          x = 'Body Mass (kg)') %>%
  gf_theme(scale_color_manual(values = my_colors)) %>%
  ggplotly(tooltip = 'text') 

cardiac_mass

htmlwidgets::saveWidget(cardiac_mass, "figures/cardiac_mass.html")

cardiac_mass_log <- gf_point(boot_median_cardiac ~ mean_mass, data = cardiac_preds,
         color = ~habitat,
         text = ~animal) %>%
  gf_errorbar(boot_lo + boot_hi ~ mean_mass, data = cardiac_preds,
              color = ~habitat) %>%
  gf_refine(scale_x_continuous(trans='log10')) %>%
  gf_refine(scale_y_continuous(trans='log10')) %>%
  gf_labs(title = 'Expected Cardiac Output\nfor the species modeled',
          y = 'Predicted Cardiac Output (mL/min)), Log10 scale',
          x = 'Body Mass (kg), Log10 scale') %>%
  gf_theme(scale_color_manual(values = my_colors)) %>%
  ggplotly(tooltip = 'text')

cardiac_mass_log

htmlwidgets::saveWidget(cardiac_mass_log, "figures/cardiac_mass_log.html")

Notes: The figure could also, of course, be non-interactive.

We can also make predictions for a “random selection of hypothetical new species of various masses” instead if you want to get a graph of lines with symmetrical error bands.

Total Ventilation

Total ventilation is the product of breathing rate (fr) and tidal volume (vt).

We will use a parametric bootstrap to make predictions of total ventilation as a function of log10(body mass) and habitat based on our random-effect models for fr and vt.

We will make sets of many (1000) predictions for each of the species that occur in both the fr and vt data sets.

# note NAMES HAVE TO BE SAME as passed to glmmTMB() when models were fitted (geesh)
mammal_fr <- readRDS('fitted-models/fr-data.RDS')
fr_mod <- readRDS('fitted-models/fr-re-model.RDS')
mammal_vt <- readRDS('fitted-models/vt-data.RDS')
vt_mod <- readRDS('fitted-models/vt-re-model.RDS')

vent_species <- intersect(pull(mammal_fr, animal), pull(mammal_vt, animal))

The two datasets have 31 species in common.

Make datasets for which to make predictions, containing these species.

vent_fr <- mammal_fr %>%
  dplyr::filter(animal %in% vent_species) %>%
  dplyr::select(animal, order, genus, species, log10.body.mass, log10fr, habitat) %>%
  mutate(across(where(is.factor), factor))

vent_vt <- mammal_vt %>%
  dplyr::filter(animal %in% vent_species) %>%
  dplyr::select(animal, order, genus, species, log10.body.mass, log10vt, habitat)

For each row of the dataset, make a set of 1000 model predictions of ventilation. These predictions will include uncertainty in the model parameters (intercept and slopes) as well as variation attributed to order/genus/species. (Residual variation is not included). The mass estimates are different in the different datasets; here we use the mass estimate from the dataset to which the model (fr or vt) was originally fitted.